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NOMENCLATURE 


A — first end-body 
B — second end-body 
C = center of mass of the tether system 
E — longitudinal stiffness of the tether 
Jc = moment of inertia of the tether system 
L — total tether length 
Lk = length of the k-th segment of the tether 
m A = mass of the first end-body 
m B — mass of the second end-body 
rftfc = embedded mass k 
R — geocentric radius 
R e = mean radius of the Earth 
R c ~ geocentric radius of the center of mass 
s = arclength along the unstretched tether 
T — tether tension 
t — time 

7 = tether elongation 
p — tether mass per unit length 
A = longitude 

T = unit vector along the tether hne 
T x — direction of the imaginary straight tether line 
Q — rotational angular rate 
( ‘ ) — differentiation with respect to time 
( 1 ) ~ differentiation with respect to the arclength 



1. FORMULATION OF THE PROBLEM 


In this study, we consider the dynamics of a spinning tether system in an 
elliptical orbit in application to the Momentum Exchange Electrodynamic Reboost 
system. Momentum exchange tether systems have been studied in a variety of ap- 
plications since Hans Moravec’s early publication [1]. It has recently been suggested 
that momentum exchange systems can be enhanced with electrodynamic reboost 
between payload transfers [2-4]. 

The Momentum Exchange Electrodynamic Reboost system (MXER) has a 
projected tether span of up to 100 km, and spins rapidly with a period of 6-7 min. 
It is placed in an orbit with a low perigee of about 400 km and a high apogee of 
about 8000 km. To capture a payload at a perigee rendezvous, within a window of a 
few seconds, the motion of the system has to be predicted with very high precision, 
having acceptable position errors on the order of 1 m. 

While this level of precision is routinely achieved today for conventional (non- 
tethered) satellites, it is much more difficult to achieve for a 100 km long flexible 
tether system. It is the goal of this study to investigate theoretical aspects of the 
dynamics and offer a practical approach to high precision dynamic modeling of a 
typical momentum exchange tether system. 


2. DYNAMIC MODEL AND EQUATIONS OF MOTION 


For the purposes of this study, we will assume that the momentum exchange 
system consists of two end-bodies A and B , modeled as point masses m A and m B , 
respectively, and a number of embedded masses m*,, k — 1, . . . , K , connected with 
tether segments of lengths Lk, as shown in Fig. 1. Mass m A can be counted as m 0 , 
and mass m B as m K +\. Segment Lk connects masses mj, and rrik-\- \ ■ 


Lq Li 


Lk 


m A m, m 2 ... m k m k+ , m B 


Fig. 1. Structure of the momentum exchange tether system. 
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All point masses and the masses of the tether segments are assumed to be 
constant. Tether mass per unit length can vary along the tether. 

Positions of the tether elements with respect to a non-rotating geocentric 
reference frame OXY Z are defined by the geocentric radius R as a function of the 
arclength s measured along the unstreched tether from A to B, and time t, 

R = R(s,t), 

Positions of the end masses and embedded masses are 

R-A = Rfl — R(SS;t)> Rfc = R(sfc,t). 


k 


A 



Fig. 2. Positions of the tether system elements. 


The tension vector of a perfectly flexible tether is tangent to the tether line 

T = — R\ 7 = l R 'l- (1) 

7 

where prime denotes differentiation with respect to the arclength s, and 7 is the 
local tether elongation. 

The tether tension can be expressed as a function of the elongation 7, elon- 
gation rate 7, temperature 0, and other factors, such as creep history, 

T = T(s, t, 7, 7, 0 ,...). ( 2 ) 


14 



Equations of motion of the tether system constitute a mix of ordinary and 
partial differential equations [5]. The motion of the end masses and embedded 
masses is described by the ordinary differential equations 

m A R^ = T a + m A g A + F a 

m B Rb — — Tb + m B g B + Fg (3) 

m k Rjt = Ta. + - T fe _ +m k g k + F k 

where dots denote differentiation with respect to time, g^, g B , and g k are the 
gravity accelerations at points A, B, and k , respectively, F A , F B , and F*. are non- 
gravitational forces acting on the end masses and embedded masses. The tether 
tension vectors are taken at the following points: T A at point A, T B at point f?, 
Tfc_ at point k of segment L k ~\ , and T*. + at point k of segment L k . 

The motion of the tether is described by the partial differential equation 

pR = T'+pg + F (4) 

where dots denote differentiation with respect to time t, and primes denote differ- 
entiation with respect to the arclength s, p is the tether mass per unit length, and 
T is the tether tension. 


3. GRAVITATIONAL FIELD MODEL 


The gravitational field is a sum of the gravitational field of the Earth and 
other celestial bodies, primarily, the Moon and the Sun. 

The geopotential U is usually represented as 


U 


I^E 

~R 


EE 

n=0 m— 0 


Re 

Ji 


{Cnm cos m\ + S nm sin mA) P nm (cos a), (5) 


where is the gravitational constant of the Earth, R E is the equatorial radius of 
the Earth, R is the geocentric radius, A is the geographical longitude eastward from 
Greenwich, a is the geocentric colatitude (the angle between the rotation axis of 
the Earth and the geocentric radius of a point, a = 0 at the North Pole), P nm are 
normalized Legendre functions, 

Pnm(cosa) = x™P™(cosa), 
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is a norm, 


= \/2n + 1, 


yt™ 1 — 


1 2(2 n + 1) {n — m)\ 
(n + m)\ 


with m = 0, 
with m > 0, 


and P™ are associated Legendre functions of degree n and order m, 


P™(*) = ( l-* 2 ) 


_ 2 \m/2 


d r 


dx m 


Pn{x ), 


P n (x) = 


d n 


2 n n! dx n 


(x> - 1)", 


P° _ 1 

-r Q — i, 

P^cosa) = cos a, 

P* (cos a) = sin a, 

^2° (cos a) = -(3 cos 2a + 1) = - (3 cos 2 a - 1), 

3 

Pi (cos a) = - sin 2a = 3 sin a cos a, 

Z 

3 

P|(cosa) = -(1 — cos 2a) = 3sin 2 a, 

Z 

P“ (cos a) = -(5 cos 3a + 3 cos a) = -(5 cos 3 a — 3 cos a), 
8 2 

3 3 

P 3 1 (cosa) = -(sin a + 5 sin 3a) = - sina(5cos 2 a — 1), 

8 2 

7~>2 / \ 1 5 , . o 

P 3 (cos a) = —(cos a — cos 3a) = 15 sin a cos a, 

15 

P 3 (cos a) = — (3 sin a — sin 3a) = 15 sin 3 a, 


The components of the gravitational acceleration g = VU are 

_ d u _ i du i du 

9R ~ dR' 9a ~ Rda' 9X ~ Psina ’ 

where <7 fl is pointing along the geocentric radius, g a is pointing southward along 
the meridian, and g\ is pointing eastward along the parallel. In the geocentric axes 
OXYZ, 

g x — ( gR cos <p + g a sin <p) cos A — g\ sin A 
9y = ( 9r cos + g a sin <p) sin A + g x cos A 
9z = 9r sin ~ g a cos <p 

where p> = tt/2 — a is the latitude. 
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To properly evaluate convergence, we should use a geopotential model of the 
maximum available degree and order. An obvious choice is the Earth Gravitv Model 
EGM96 [6], of degree and order 360. 7 

According to this model, fi E = 398600.4415 km 3 /sec 2 , R E = 6378.1363 km, 
and the normalized coefficients C nm and S nm are as follows: 


C'o.o = 1, 

C'i.o = 0, 

^ 1,1 = 0, 

C 2 ,o = -0.484165371736 • 10~ 3 , 
C 2 ,i = -0.186987635955 • 10“ 9 , 
C 2 ,2 = 0.243914352398 • 10~ 5 , 
C 3 , o = 0.957254173792 • 10 -6 , 
C 3 , i = 0.202998882184 • 10" 5 , 
C 3 , 2 = 0.904627768605 • 10~ 6 , 
C 3) 3 = 0.721072657057 • 10" 6 , 


50.0 — not used 
-S'ljO — not used 
§1,1 = 0 , 

S 2t o — not used 

5 2 .1 = 0.119528012031 • 10 -8 , 
S 2l2 = - 0.140016683654 • 10 ~ 5 , 
S 3< o — not used 

§z , i = 0.248513158716 • 10 -6 , 

5 3 . 2 = - 0.619025944205 • 10 “ 6 , 

5 3 . 3 = 0.141435626958 • 10 ~ 5 , 


<^ 360,359 = 0.183971631467 • 10 10 , S 36 0,359 = -0.310123632209 • 10 _1 °, 

^360,360 = -0.447516389678 • 10 -24 , 5 36 o,36o = -0.830224945525 • 10 -1 ° 

Gravitational perturbations caused by the Earth tides can be described by 
time dependent variations in the geopotential expansion coefficients [7], 

Gravitational perturbations caused by a celestial body of mass M p located at 
R p with respect to the geocentric axes OXY Z are described as 


gp — GM p 


/ Rp — R 

VIRp-RI* 



( 6 ) 


where G is the universal gravity constant, and R v = |R p |. 

For most practical purposes and for all celestial bodies, except the Moon, the 
linear approximation 

GM P . , 

gp — [3e p (ep,R) — R], (7) 

where e p = R P /R P is a unit vector of the direction to the celestial body, should be 
adequate [7]. 
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4. STATIONARY ROTATION 


If we disregard non-gravitational forces F in equations ( 3 )-( 4 ) and assume 
that the tether is inextensible (7 = |R'| = 1 ) and gravitational acceleration does 
not vary along the tether (g = g c ), then equations ( 3 )-( 4 ) will have a stationary 
solution 

T = T(s)t(2), R — Rc(0 + (<$ “ 6 C ) t(<), (8) 

where R c(t) and t(£) are determined from 


Rc = gc, T = —Q 2 T. 


(9) 


Here, R c is the radius- vector of the center of mass, moving as a point mass, g c is 
the gravity acceleration at the center of mass, T is a unit vector along the tether line 
rotating at a constant angular velocity Q in a fixed plane normal to O, Q = |Q|, 
s c is the arclength corresponding to the center of mass 


s c 


— — (^rnA s A + r fnB s B + ^ ] mk s k + J p s ds^j , 


( 10 ) 


and M is the total mass of the tether system. 
In the stationary motion ( 8 )— (9) , we have 


R = R c + (s - s c ) T = gc - (-s - s c ) f? 2 T. 


After substituting this relation into equations (3)-(4), we arrive at a boundary 
problem for the equilibrium tension, 

T' = —pfi 2 [s — Sc), 

T a = —Tn A 0 2 (sA — Sc), 

Tb = ttib&{sb Sc), 

Tk+ = Tk- — mfcf? 2 (sfc — s c ), 

where the last equation defines tension increments associated with the embedded 
masses m^. 

The equilibrium tension is obtained by solving the first equation of (11) with 
the initial condition T = T A given by the second equation and incrementing the 
tension at the points Sk, as described by the last equation of (11). The boundary 


18 



condition T — T B given by third equation is always satisfied because of the definition 
of the value s c (10) and the structure of equations (11). 

In general, the equilibrium tension increases from the value T A at the end A 
to its maximum value Tc at the center of mass C , and then drops to the value Tg 
at the other end, as shown in Fig. 3. 

To minimize mass, the tether should be tapered, with the linear density vary- 
ing along the tether, p — p(s). Ideally, the tether should be equally stressed/strained 
at all points, 


T 

E 


= 8 




( 12 ) 


where E is the longitudinal tether stiffness, and is the maximum allowed strain. 
Condition (12) can be rewritten as 


T = 


P V E ? 



(13) 


where v E is the longitudinal wave velocity in the tether. Now, we can express the 
linear density as a function of tension, p = Tv^ 2 S * 1 , and substitute this relation 
into the first equation of the boundary problem (11) to produce a boundary problem 
for an ideal equally stressed/strained tether. 


t (km 



Fig. 3. Tension profile of a tether system with (a) and without (b) payload. 




Ill practice, the tether can be made of a number of uniform segments, approx- 
imating the desired taper profile. 

Typical tension profiles are shown in Fig. 3 for a 90 km tether system with 
m A - 250 kg, m B - 11000 kg, and 8 embedded masses m k = 200 kg placed 10 km 
apart. The tethers are made of Zylon and are uniform on each segment with a 
maximum strain of S„ = 0.01. The system spins at 12 = 0.8 deg/sec. In case (a), 
the end mass A carries a 2500 kg payload, and there is no payload in case (b). 


5. SMALL OSCILLATIONS AND EIGENFORMS 


To study small oscillations of the tether system about the stationary rotation, 
we introduce a rotating reference frame C£r](, with the origin at the center of mass, 
C£ aligned with the tether line, and axis G aligned with the angular rate vector 
n, so that 

t= (1,0,0), Q = (0, 0, Q). 

With respect to the rotating axes C£rj(, the stationary motion (8)-(9) is viewed as 
a relative equilibrium. 

Equations of small oscillations of the tether system about the relative equi- 
librium are derived from (3)-(4), 

pD(8 r) = 6T', 
m A D(8r A ) = 6T a , 
m B D(8 t b ) = -8T b , 
m k D(8v k ) = 6T k+ - 

where 8r and 8 T are deviations from the relative equilibrium, and D denotes the 
linear expression 

D(8r) = 8r + 2D. x 8r + Q x (D x 8r), 
in which derivatives are calculated with respect to the rotating axes C£r](. 

Small oscillations (14) have two independent components. One is normal to 
the tether line in the plane of the stationary rotation, 

£r = (0, 77, 0), 8T = (0, Tr, ', 0), 

and the other one is normal to the rotation plane 

8v = (0, 0, C), ST = (0, 0, T('). 


There is no longitudinal component in the linear approximation for an inextensible 
tether. 



Equations of small in-plane oscillations take the form 


p{r,-Q 2 r ] )^{T7 1 , )\ 
m A (rj A - f2 2 r] A ) = (Tt]') A , 
m B (vb - (2 2 Vb) = -(Trj') B , 
m k (rjk ~ fi 2 Vk) = (Trj') k + - (Tr]') k -, 


( 15 ) 


After substituting rj = U n (s) cos(Q n t) into (15), we arrive at the following eigen 
value problem, 

(Tu'j = - P (nl + tf)u n , 

{Tu' n ) A = -m A {n 2 n + n 2 )u nA , 

(TU’ n )a = m B + 12 2 ) U nB , 

(TV n ) k+ = ( TU’ n ) k - - m k (ill + Q 2 ) U nk , 

Small out-of-plane oscillations are described by 


(16) 


p< = (T(’y, 

m A Ca = (T (') A , 

t^bCb = —{T(')b, 
m k ( k =(T(') k+ -(T(') k -. 


(17) 


After substituting ( = U n (s ) cos(fl n t) into (17), we derive the following eigenvalue 
problem, 

(tu'J = - P nlu n , 

(TUUa = -m a alu nA , 

(TU’ n ) B = m B alU nB , 

{TU' n ) k+ = (TU' n ) k --m k nl U nk . 

The eigenvalue problems (16) and (18) are very similar. They have the same 
eigenfunctions U n (s ), and their eigenfrequencies are bound by a simple relation, 


n 2 n< = nL + n 2 . 


717 ? 


(19) 


Analyzing the boundary problem (11) for the equilibrium tension and the 
eigenvalue problems (16), (18), we note that the eigenforms do not depend on the 
angular rate !?. Indeed, we can introduce a normalized equilibrium tension 


p M = 


t W 
n- 


( 20 ) 


21 



and rewrite equations (11) as 


P' = -p(s - s c ), 

Pa — ~ rr >-A{ s A — •Sc)) 

P B = m B (s B - s c ), 

Pk+ = Pk- ~ rn k (s k ~ s c ), 

while equations (16) and (18) can be reduced to 

(PU'J = -pftU n , 

( pu L)a = -m A plU nA , 

{PU' n ) B = rn B (3l U nB , 

( PU'n)k+ = ( PU'Jk - - m k Pi Unk, 


(21) 


( 22 ) 


where P n are dimensionless eigenvalues. The resulting system of equations (21)- 
(22) and its solutions do not depend on the angular rate Q. 

The eigenfrequencies Q nv and Q n ^ can be expressed through the dimensionless 
eigenvalues ft n as 


n nv = y/Pi-i n, n n( =p n n. ( 23 ) 

The eigenvalue problem (21)-(22) has two trivial solutions, 

Po = 0, U 0 = 1, (24) 

and 

P\ = 1) Ui = s - s c . (25) 

The first solution simply reduces the right and the left parts of (22) to zero, while 
the second solution reduces (22) to (21). Dynamically, solution (24) corresponds to 
perturbations of the motion of the center of mass, while solution (25) corresponds 
to variations of the angular rate and direction of the axis of the stationary rotation. 

The trivial solutions are followed by an infinite series of solutions {P n , U n }, 
n = 2, 3, ... , which constitute an orthogonal basis with the orthogonality condition 

(Ui, Uj) = 0, i j, (26) 


where 


(Vi, Vi) 


— m A(UiUj)A + TnsiUiU j)s + 


J2 m ^iUj) k 

k 



p(UiUj) ds. 


(27) 
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23 










( 32 ) 


Also, the following relation for the eigenvalues can be derived from (22) 

0l = ^~J*P(K) 2 ds = ci. 

Some eigenforms are shown in Fig. 4-5 for the same parameters as in case (b) 
of Fig. 3, namely, L = 90 km, m A = 250 kg, m B = 11000 kg, m k = 200 kg, k = 
1,. . . ,8, connected with 10 km long uniform tethers of the following masses: 318, 
365, 406, 438, 461, 473, 475, 473, 460 kg. The corresponding nontrivial eigenvalues 
(re = 2, ... , 19) are given below 

P 2 = 2.136, /3 a = 10.497, /? 14 = 20.182, 

Pz = 3.339, p 9 = 11.355, p 1B = 20.573, 

Pa = 4.590, p 10 = 12.414, p 16 = 22.015, 

Ps = 5.912, p^ = 15.549, p 17 = 23.714, 

Pe = 7.328, p 12 = 18.020, p ls = 25.550, 

P 7 = 8.853, p 13 = 19.320, p 19 = 27.529. 

Generally, the n-th form has re nodes. Note that while the lower modes (Fig. 4) 
show smooth curves of average collective behavior, the higher modes (Fig. 5) exhibit 
a more peculiar interplay between the tether and embedded mass dynamics. 

6. DECOMPOSITION OF MOTION 


Using eigenforms obtained from the solution of the boundary problem (21) 
and the eigenvalue problem (22) for a simplified case of a stationary rotation, we can 
represent solutions of the general equations (3)-(4) in the form of a series expansion 

OO 

R ( s > <) = X! **"(*) U M> (33) 

n—0 

where the term re = 0 (trivial form U 0 = 1) represents, in essence, the orbital motion 
of the center of mass, q 0 = Rc(t), and the term re = 1 (trivial form U\ = s — s c ) 
describes a quasi-rigid rotation of an imaginary straight tether, q x = T 1 (t), where 
is pointing along the imaginary straight tether line, so that 

OO 

R(s, t) = Rc(0 + (5 — s c ) Tj (t) + ^2 q n(t) U n (s). (34) 

71=2 
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The higher order terms n > 2 represent tether oscillations about the imaginary 
straight line drawn through the center of mass C along vector t x . Kinematically, 
they do not contribute to the displacement of the center of mass and the quasi-rigid 
rotation because of the orthogonality conditions (26). 

Note that the term with T x also describes a uniform elongation of the tether 
system. If we disregard the higher order terms n > 2, then the tether elongation 
7 = |R'| = | T i| = 7 i will be the same for all tether elements. 

Expression (33) yields similar kinematic and dynamic relations, 

00 oo 

R( 5 , t) = J2 qn(<) U n (s), R (M) = Y, q n(<) U n ( S ), (35) 

n=0 n=0 


or in terms of (34), 


R(s, t) — R c (t) + (-5 s c ) Ti(t) + 'Y, Qn{t)U n (s), 

71=2 

oo 

R(a, i) = R c (i) -f (if — Sc)Ti(() + <in(t)U n (s), 


(36) 


71=2 


where the derivatives are calculated with respect to the non-rotating geocentric 
frame OXYZ. 

We can now substitute the second equation of (35) into the general equations 
(3) (4), multiply by U n , n — 0, 1, 2, . . . , and sum over all tether elements and point 
masses. Applying the orthogonality conditions (26), we find that 


Qn, — 


Wn\ 


(Qra 4" + ®n)) 


n = 0,1,2,..., 


(37) 


ds 


where Q n is reduced to an integral of tension, 

f B 

Qn = T A U nA -T B U nB + Y( T k+-T k -)U nk + / T'Z7 n 

k 

= - f B T V n ds, 

J A 

G n is expressed through the inner product (27), 

G» = (g, U n ) 

f B 

= m A (gU n ) A + m B (gU n ) B + Y m k(gU n ) k + / p(gU n )ds , 

L j A 


(38) 


(39) 



and <D n is a generalized sum of non-gravitational forces, 

f B 

^aUtia + FbUub + ^ ] Fj, U n k + / F U n ds. (40) 

k J A 


The system of ordinary differential equations (37) is equivalent to the original 
system of mixed partial and ordinary differential equations ( 3 )-( 4 ), because it is 
derived through the general transformation of variables ( 33 ) without any simplifying 
assumptions. 

Using notations (34), and taking into account that ||f7 0 || = M and Q 0 vanishes 
because U 0 = 0 , the motion of the center of mass is described by 

R-c = — (G 0 + <D 0 ). (41) 


Taking into account that ||17i|| — J c and U[ — 1 , the quasi-rigid rotation and 
uniform elongation are described by 


Ti = 


(Qi + Gi + ®i), 


Q 


i — 


/: 


T ds. 


According to (34), the tangent to the tether line is defined as 

OO 

R' = t 1 + ^t, Ar= '£q n U' n . 

71=2 

Assuming that |utJ <£ T, , , the elongation can be approximated as 


(42) 


(43) 


T= l R 'l = V / Ki| 2 + 2(T 1 ,aT) + |ziT| 2 « 7l +(e l ,4T) = 7 1 + ^T qen U' n , (44) 

71=2 

where 

Tl K|, 7 9en, = (®lj*fn)' (4b) 

The tether tension depends not only on the elongation, but also on temper- 
ature, internal friction, creep history, and other factors. The following expression 
can be adopted within a certain range of conditions, 


T = T 0 + E(j - 7 *) + Dj, (46) 

where T 0 is a reference tension, E is the tether longitudinal stiffness, D is the 
damping coefficient, 

7 * = l0 + e(t-t 0 ) + a(G-e 0 ), ( 47 ) 
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7 o is the elongation under the reference tension T 0 at time t 0 and a reference tem- 
perature 6 > 0 , e is the creep rate, a is the thermal expansion coefficient, and 0 is the 
current temperature of the tether element. Other terms can be added as needed. 
Substituting (44) into (46), we derive 

oo 

T — Ti q en + D q en ) U n , T t = T 0 + E (ji — 7 *) + (48) 

n~ 2 

and then calculate 


T 

7 




9en + D q en 


UL 


(49) 


Retaining only linear terms in q n , the integral of tension Q„ can now be 
represented as 


Qr 


T T / 00 

= ~Ja 7 R,U n is = -J A “( T > + X>tfi)tC<fa 

k~2 

fB 00 

* -J ( t ^ + T, c ‘‘U' kjKds, 


(50) 


where ej — T x / 71 is a unit vector along the imaginary straight tether line and 


T 

c n = Ee^qen + De 1 q en + ~ (q„ - e l9en ), 

7i 


?en = (e 1 ,q n ). (51) 


Expression (50) can be significantly simplified in the ideal case of a perfectly 
tapered tether. Let us define the tether reference state T 0 (s) as a stationary rotation 
at an angular rate f2 0 free of external torques, as described in Section 4 , with the 
exception that the tether is now elastic. Let us assume that the tether is perfectly 
tapered so that all tether elements have the same elongation 7 0 under the stationary 
tension T 0 at a given temperature 0 O . Then, similarly to ( 11 ), the equilibrium 
tension profile can be determined from the boundary problem 

T' = -pfi 2 (s - s c ) 70 , 

Ta = —m A Q (5^4 — Sc) Jo, 

T B = m B fi 2 (s B - s c ) 70 , ^ 

Tk-\- Tfc— (<S/fc Sc) 7o , 

and similarly to ( 20 ), the equilibrium tension can be expressed as 

To(s) = f2 0 7o P($), (53) 
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where P(s) is the solution of the boundary problem ( 21 ). 

Within the linear elasticity region, the longitudinal stiffness of the tether can 
be related to its tension as 

B(s) = ^\=k. flj P(s), = * (54) 

With a strain on the order of 7o - 1 ~ 0.01, the coefficient k E ~ 100. 

The internal damping coefficient is usually considered directly proportional to 
the longitudinal stiffness and inversely proportional to the frequency of oscillations 
f? n , which gives 


D(s,f2 n ) 


n n 


v n n 0 p(s), 


Vn 


k E rj 

^T’ 


(55) 


where tj is a loss factor, on the order of 0.1 for braided tethers. 

After these substitutions, the first term under the integral in formula (50) 
takes the form 


T iei = QlP{s )[ l0 + k E { lt -7*) + i/i7i] ei = Q 2 0 P{s)u 1 x u 

where 

— [Tq+^Ct 1 -TO + ^iTi], (56) 

and the expansion coefficients become 

Cn = f2 0 P(s) [& E e 1 q en + u n exg en + (q n — e^n)] . 

If we assume that 7 * defined by (47) does not vary along the tether, i.e., the 
creep rate, thermal expansion coefficient, and the temperature are the same for all 
tether elements, then we can apply the orthogonality condition (28) along with the 
eigenvalue relation (32) to derive 

Qi = -Jo^u 1 x 1 , ( 57 ) 

for n = 1 and 

Q™ ~ — ll^nll [kE Gi9en + V n e 1 q en + U x (q n — e x ^erj.)] (58) 

for n — 2 , 3 ,... 

Now, we can rewrite the equation of quasi-rigid rotation (42) as 



and the equations of tether oscillations (37), n = 2 , 3 ,..., 


as 


q n +/?n^O [ fc B e l?en + ^nej^en +«l (q„ - e^en)] = -J_ (Q n + <D n ). (60) 


7 . GRAVITATIONAL FORCES 


The gravitational terms in the equations of motion (41), (59) and (60) are de- 
fined by (39). The variation of the main component of the gravitational acceleration 
g along the tether can be represented in the linear approximation in displacement 
from the center of mass as 


g = 


He R 



He 


[3e* (e B ,r) - r], 


(61) 


where He is the gravity constant of the Earth, e R = R C /R C is a unit vector along 
the geocentric radius of the center of mass R c , and 


r( 5 , t) = (s- s c ) t x ( t ) + q ra (t) U n (s ) (62) 

n = 2 

is the position vector relative to the center of mass. 

Substituting expressions (61)-(62) into (39), and applying the orthogonality 
conditions (27) and the norm definition (29), we obtain 


G 0 

G x 


- M 



J o [^(en.T,) — T t ], 

\\Vn\\^ [3 ejt ( eit ,q„)-q„], 


n ~~ 2,3,..., 


(63) 


where M is the total mass of the system (30), and J c is the moment of inertia of 
an imaginary system with an unstretched straight tether about the center of mass 
defined by (31). 

Note that in the linear approximation, the higher modes of oscillations do 
not affect the motion of the center of mass and the quasi-rigid rotation about the 
center of mass, and there is no practical need to include higher order terms in q n as 
long as typical amplitudes of tether oscillations in a momentum exchange system 
are relatively small. However, to achieve the required precision, it is necessary to 
include higher order terms in (s — s^)^. 
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The variation of the main (Newtonian) term of the gravitational acceleration 
along the line ( s - s c ) x 1 drawn through the center of mass is given by 

g — ^ E ^ _ _l l E [Rc + (s — ^c)^l] 

R 3 |Rc + (5 — S c ) Ti | 3 

where 

* = -(e«,e 1 ), e x = — , 

h 

As known in the theory of Legendre functions, 


(1 - 2are + e 2 ) 1 / 2 X P n(x) 

7 n-0 

where P n (x) are Legendre polynomials. Differentiating this equation with respect 
to x and dividing by e, we find that 


— _ f J ' E e n T ee i , 

R 2 C (1 - 2xe + e 2 ) 3 / 2 ^ 


( s ~ fo) 7i 


(l - 2x£ + £ 2 ) 3 / 2 X^+lC*)^ 

7 n=0 


where -P n are the first derivatives of the Legendre polynomials. 
After substituting this expression in (64), we obtain 


g ~~W X i e * P n+i( x ) + e i K(x)] 


Rl 




n=0 


(* - g c) 7i 


(65) 


This series expansion converges rapidly because ( s-s c )/R c is typically small, 
on the order of 0.01. To determine a minimum number of terms required for 
practical computation, consider a simple inequality 


n 



or 


> In (n/S) 

~ In {Rc/LY 


(66) 


where the left part estimates the contribution of the n-th term, L is the total tether 
length, and 6 is a maximum acceptable relative error. Generally, all terms N, 
where N is the minimum number satisfying (66), can be dropped. For a typical 
momentum exchange system, it is sufficient to retain only 5 terms at the perigee, 
and 4 terms at the apogee to achieve 10-digit precision ( 8 lO” 10 ). 

Expansion (65) can be also presented as 


g = Xsc(*-sc) n , 

n— 0 


K » = if!Zs 

c n\ d s n 


( 67 ) 
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where the derivatives with respect to s are calculated along the line (s - 3c ) Tl 
drawn through the center of mass. For the Newtonian term, according to (65), 

= [ e *-Pn+l(*) + ei P' n {x)\ (68) 

Note that expansion (67) can be treated more broadly to include all terms of 
the standard representation of the Earth gravitational field, as well as the gravita- 
tional fields of other celestial bodies. 

Using (67), the generalized gravitational forces G n (39) can be represented as 


— (( s ~ s c) k , U n ). (69) 


Applying the orthogonality conditions (27), norm definition (29), and trivial 
form properties (24)-(25), (30)-(31), we find that 



where I n are the high order moments, 


/„* = h, 

If =Ik + l, 


(70) 


In 


m A(sA s c ) n + m B (s B — 5 C ) n -(- 'y ^ mfc (sfc 

k 



s c ) n ds. (71) 


We note also that 


r 

n 


In = 0 , 


n = 2,3,. 


(72) 


The practical benefit of using expansion (69) is that the quantities can be 
precomputed to avoid integration of the gravitational force along the tether on each 
step of solution propagation. 

Another practical question is to estimate how many terms of the gravitational 
field expansion (5) must be retained to achieve the required accuracy. The problem 
is in the slow convergence of the series expansion at low altitudes when the ratio 
R E / R is close to 1. 


On average, according to the Kaula rule [7], the normalized geopotential co- 
efficients C nrn and S nrn decrease in inverse proportion to the second power of their 
degree, 


Cnm ? $r 


10 - 


n * 


(73) 
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o 1 2 t (h) 3 


Fig. 6. Minimum degree of the geopotential expansion for MXER. 


Taking this empirical rule into account, the following inequality roughly esti- 
mates the contribution of the higher degree terms to the acceleration of the center 
of mass of a momentum exchange tether system, 


10“ 5 

n 



^ In (nS • 10 s ) 

U ~ ]n(R c /R E ) ’ 


(74) 


Numerical experiments show that a value of 6 ~ 10 -9 can be assumed for 
MXER to achieve better than 1 m prediction accuracy over one orbit. The minimum 
degree N that satisfies (74) is shown in Fig. 6 as a function of the geocentric radius 
and flight time. We see that a lower degree model is adequate for the most part of 
the orbit, but a model of a much higher degree is needed during perigee passages. 
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8. SIMULATION APPROACH 


Now, we can gather all equations, and apply formulas for the gravitational 
forces derived m the previous section. The equation of motion of the center of : 

(41) takes the form 


mass 


®o 

Af ’ 


k- 2 

the equation of the quasi-rigid rotation (59) can be rewritten 

1 °° 

+ fil Ul x, = gl + -1 + p. 


(75) 


as 


C 1^2 J ° 


(76) 


and the equations of tether oscillation (37), n = 2, 3, . . . , can be reduced to 

qn+ [k E ®1 9en + (q n — e 1 g' en )] = 

g(, 


(77) 


k~ 2 


As noted earlier, very few terms g£ are needed to achieve the desired precision, 
and quantities can be precomputed. 

Equations (75)-(77) cleanly separate three dynamicahy different kinds of mo- 
tions of the momentum exchange tether system. They are compact and flexible 
in sense that the number of the gravitational terms and the number of the tether 
oscillation modes can be selected depending on the desired accuracy. This opens a 
way to building a very computationally effective simulator. 

Preliminary testing for a typical momentum exchange system showed that 
equations (75)-(77) can be integrated with a time step of about 4 seconds of orbital 
time, and it takes only a few seconds on a standard issue PC to propagate a solution 
over one orbit with 9-10 digit precision. More studies are needed to understand 
possible limitations of this approach. 
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